function calcCellParamsComWind
% calcCellParamsComWind
%
% Basically a temporary kluge to calc cell params for traces 
% acq'd with older version of ephus. 
%
% To use, click on the ephys trace. 
%
% NB: assumes the following pulse param's were used:
%                 sampleRate: 10000
%                  amplitude: -5
%                     offset: 0
%     squarePulseTrainNumber: 1
%        squarePulseTrainISI: 0.2000
%      squarePulseTrainWidth: 0.1000
%      squarePulseTrainDelay: 0.1000
%
% Note that Ri and Rm are redundant.
%
% gs 2008 04
% ----------------------------------

% trace data
traw = get(gco, 'YData')';

% stimulus parameters
sr = 10000;
V = -5;
numsteps = 1;
isi = .1;
delay = .1;
dur = .1;

% calculate baseline and baseline-subtract the whole trace
snippetAtStart = traw(1 : delay * sr - 1);
snippetAtEnd = traw((end - round(0.02 * sr)) : end); % last 20 msec
baseline = mean([snippetAtStart; snippetAtEnd]);
tbs = traw - baseline;

% calculate an average response
for n = 1 : numsteps
    startInd = round((delay + (n-1) * isi) * sr);
    endInd = startInd + round(dur * sr);
    tracelet(:,n) = tbs(startInd : endInd);
end
tmeanbs = mean(tracelet,2);

% make a corresponding time trace
x = 0:length(tmeanbs)-1;
x = x / sr;

% find peak, within first 10 msec
[Ipeak, peakX] = min(tmeanbs(1:round(0.01*sr)));
% assumes pulse is hyperpolarizing; response is downward

% steady state: average over last 20 msec of response
steadystate = mean(tmeanbs(end-round(0.02*sr):end));

% calculate Rseries and Rinput
Rs = 1000 * V / Ipeak;
Ri = 1000 * V / steadystate - Rs;

% steadystate-subtracted mean tracelet:
tmeanss = tmeanbs - steadystate;

% calculate Cm, based on time for signal to fall to 1/e of peak
efold = find(tmeanss(peakX:end) >= (Ipeak-steadystate)/exp(1));
dropInd = efold(1);
tau = x(dropInd);
Cm = 10^6 * (Rs + Ri) * tau / (Rs * Ri);

% graphics
figure('Color', 'w'); 
subplot(2,1,1)
plot(x,tmeanbs, 'b-')

subplot(2,1,2)
plot(x(peakX:end), tmeanss(peakX:end));
hold on
plot(tau, tmeanss(dropInd),'ro');

% command line display
disp(' ');
disp(['map.Rs = ' num2str(round(Rs)) ';']);
disp(['map.Ri = ' num2str(round(Ri)) ';']);
disp(['map.Rm = ' num2str(round(Ri)) ';']);
disp(['map.Cm = ' num2str(round(Cm)) ';']);
disp(['map.tau = ' num2str(tau) ';']);
